/*
 * To change this template, choose Tools | Templates
 * and open the template in the editor.
 */

package MeshlessShapeFunctions;

import HashMatrix2.HashMatrix2;
import HashMatrix2.LUFactorization;
import Integration.IFunction2D;

import java.awt.Dimension;
import java.util.ArrayList;

import javax.swing.JFrame;

import org.jzy3d.chart.Chart;
import org.jzy3d.chart.controllers.mouse.ChartMouseController;
import org.jzy3d.colors.Color;
import org.jzy3d.colors.ColorMapper;
import org.jzy3d.colors.colormaps.ColorMapRainbow;
import org.jzy3d.factories.JzyFactories;
import org.jzy3d.maths.Range;
import org.jzy3d.plot3d.builder.Builder;
import org.jzy3d.plot3d.builder.Mapper;
import org.jzy3d.plot3d.builder.concrete.OrthonormalGrid;
import org.jzy3d.plot3d.primitives.Shape;
import org.jzy3d.plot3d.primitives.axes.AxeFactory;
import org.jzy3d.plot3d.primitives.axes.ContourAxeBox;
import org.jzy3d.plot3d.primitives.axes.IAxe;
import org.jzy3d.plot3d.rendering.canvas.Quality;
import org.jzy3d.plot3d.rendering.legends.colorbars.ColorbarLegend;

/**
 *
 * @author mark_o
 */
public class JMeshlessInterpolationFunction2D implements IFunction2D{
    protected int PolinomialDegree;
    protected double xii;
    protected double yii;
    protected double rho;
    protected IKernelFunction func;
    protected HashMatrix2 A;
    protected double[] b;
    protected ArrayList<Particle> NodeList;
    double EPS = 1e-5;

    public JMeshlessInterpolationFunction2D(int PolinomialDegree, double xii, double yii, double rho, IKernelFunction func) {
        this.PolinomialDegree = PolinomialDegree;
        this.xii = xii;
        this.yii = yii;
        this.rho = rho;
        this.func = func;
    }

    public void setNodeList(ArrayList<Particle> NodeList) {
        this.NodeList = NodeList;
    }

    public double getRho() {
        return rho;
    }

    public double getXii() {
        return xii;
    }

    public double getYii() {
        return yii;
    }

    public ArrayList<Double> getBaseVector(double x,double y,double xi, double yi){
        ArrayList<Double> L = new ArrayList<Double>();
        for (int i = 0; i <= this.PolinomialDegree; i++) {
            for (int j = 0; i+j <= this.PolinomialDegree; j++) {
                L.add(Math.pow((x-xi)/rho, i)*Math.pow((y-yi)/rho, j));
            }
        }

        return L;
    }

    public ArrayList<Double> getBaseVectorDerivatyX(double x,double y,double xi, double yi){
        ArrayList<Double> L = new ArrayList<Double>();
        for (int i = 0; i <= this.PolinomialDegree; i++) {
            for (int j = 0; i+j <= this.PolinomialDegree; j++) {
                L.add(Math.pow((x-xi)/rho, i)*Math.pow((y-yi)/rho, j));
            }
        }

        return L;
    }

        public ArrayList<Double> getBaseVectorDerivatyY(double x,double y,double xi, double yi){
        ArrayList<Double> L = new ArrayList<Double>();
        for (int i = 0; i <= this.PolinomialDegree; i++) {
            for (int j = 0; i+j <= this.PolinomialDegree; j++) {
                L.add(Math.pow((x-xi)/rho, i)*Math.pow((y-yi)/rho, j));
            }
        }

        return L;
    }

    public double calculate(double x,double y){
        ArrayList<Double> BaseVector = getBaseVector(x,y, xii,yii);
        int n = BaseVector.size();
        b = new double[n];
        A = new HashMatrix2(HashMatrix2.ARRAY_TYPE, n);
        for (int i = 0; i < NodeList.size(); i++) {
            Particle particle = NodeList.get(i);
            calculateAlfa(x,y,particle.x,particle.y);
        }
        b[0]=1.0;

        double returnValue = 0.0;
            LUFactorization lu = new LUFactorization(A);
            double[] res = lu.calculate(b);


            for (int i = 0; i < b.length; i++) {
                returnValue+=BaseVector.get(i)*res[i]*func.evaluateSque((x-xii)/rho,(y-yii)/rho);
            }
        return returnValue;
    }


    private void calculateAlfa(double x,double y, double xi, double yi){
        ArrayList<Double> BaseVector = getBaseVector(x, y, xi, yi);
        int n = BaseVector.size();
        for (int i = 0; i < n; i++) {
            double val = BaseVector.get(i)*func.evaluateSque((x-xi)/rho,(y-yi)/rho);
            for (int j = 0; j < n; j++) {
                A.addValue(i, j, val*BaseVector.get(j));
            }
        }
    }

    public double calculateDiffuseDerivatyX(double x,double y){
        ArrayList<Double> BaseVector = getBaseVector(x,y, xii,yii);
        int n = BaseVector.size();
        b = new double[n];
        A = new HashMatrix2(HashMatrix2.ARRAY_TYPE, n);
        for (int i = 0; i < NodeList.size(); i++) {
            Particle particle = NodeList.get(i);
            calculateAlfa(x,y,particle.x,particle.y);
        }

        boolean key=false;
        for (int i = 0; i < b.length; i++) {
            b[i] = BaseVector.get(i)*func.evaluateSque((x-xii)/rho,(y-yii)/rho);
            if(b[i]!=0.0){key = true;}
        }

        double returnValue = 0.0;
        if(key){
            LUFactorization lu = new LUFactorization(A);
            double[] res = lu.calculate(b);

            //ArrayList<Double> BaseVector1 = getBaseVectorDerivaty(x, xii);
            for (int i = 0; i < b.length; i++) {
                returnValue+=BaseVector.get(i)*res[i]*func.evaluatePrimeSqueX( (x-xii)/rho, (y-yii)/rho);
            }


            double[] MxAlpha = new double[n];
            double tmp = 0.0;
                for (int j = 0; j < BaseVector.size(); j++) {
                    tmp += BaseVector.get(j)*res[j];
                }
            for (int i = 0; i < MxAlpha.length; i++) {
                MxAlpha[i] = tmp*BaseVector.get(i)*func.evaluatePrimeSqueX((x-xii)/rho, (y-yii)/rho);
            }

            key=false;
            for (int i = 0; i < b.length; i++) {
                b[i] -= MxAlpha[i];
                if(b[i]!=0.0){key = true;}
            }
            if(key){
                double[] res2 = lu.calculate(b);
                for (int i = 0; i < b.length; i++) {
                   returnValue+=BaseVector.get(i)*res2[i]*func.evaluateSque((x-xii)/rho, (y-yii)/rho);
                }
            }


        }
        return returnValue;
    }

    public double calculateDiffuseDerivatyY(double x,double y){
        ArrayList<Double> BaseVector = getBaseVector(x,y, xii,yii);
        int n = BaseVector.size();
        b = new double[n];
        A = new HashMatrix2(HashMatrix2.ARRAY_TYPE, n);
        for (int i = 0; i < NodeList.size(); i++) {
            Particle particle = NodeList.get(i);
            calculateAlfa(x,y,particle.x,particle.y);
        }

        boolean key=false;
        for (int i = 0; i < b.length; i++) {
            b[i] = BaseVector.get(i)*func.evaluateSque((x-xii)/rho,(y-yii)/rho);
            if(b[i]!=0.0){key = true;}
        }

        double returnValue = 0.0;
        if(key){
            LUFactorization lu = new LUFactorization(A);
            double[] res = lu.calculate(b);

            //ArrayList<Double> BaseVector1 = getBaseVectorDerivaty(x, xii);
            for (int i = 0; i < b.length; i++) {
                returnValue+=BaseVector.get(i)*res[i]*func.evaluatePrimeSqueY( (x-xii)/rho, (y-yii)/rho);
            }


            double[] MxAlpha = new double[n];
            double tmp = 0.0;
                for (int j = 0; j < BaseVector.size(); j++) {
                    tmp += BaseVector.get(j)*res[j];
                }
            for (int i = 0; i < MxAlpha.length; i++) {
                MxAlpha[i] = tmp*BaseVector.get(i)*func.evaluatePrimeSqueY((x-xii)/rho, (y-yii)/rho);
            }

            key=false;
            for (int i = 0; i < b.length; i++) {
                b[i] -= MxAlpha[i];
                if(b[i]!=0.0){key = true;}
            }
            if(key){
                double[] res2 = lu.calculate(b);
                for (int i = 0; i < b.length; i++) {
                   returnValue+=BaseVector.get(i)*res2[i]*func.evaluateSque((x-xii)/rho, (y-yii)/rho);
                }
            }


        }
        return returnValue;
    }

    public double calculateDerivatyX(double x, double y) {
    	ArrayList<Double> BaseVector = getBaseVector(x,y, xii,yii);
        int n = BaseVector.size();
        b = new double[n];
        A = new HashMatrix2(HashMatrix2.ARRAY_TYPE, n);
        for (int i = 0; i < NodeList.size(); i++) {
            Particle particle = NodeList.get(i);
            calculateAlfa(x,y,particle.x,particle.y);
        }
        b[0]=1.0;
        boolean key = true;
        double returnValue = 0.0;
        if(key){
            LUFactorization lu = new LUFactorization(A);
            double[] res = lu.calculate(b);
            
            double val =  func.evaluatePrimeSqueX((x-xii)/rho,(y-yii)/rho);//  func.evaluateDerivaty((x-xii)/rho);
            key=false;            
            double[] tmp = new double[b.length];
        	for (int j = 0; j < this.NodeList.size(); j++) { // по вузлах
        		Particle xj = NodeList.get(j);
        			val = func.evaluateDerivaty((x-xj.x)/xj.rho);
        			BaseVector = this.getBaseVector(xj.x,xj.y, xj.x,xj.y);//getBaseVector(xj.x, xj.x);
        			double tmpVal = 0.0;
        			for (int k = 0; k < res.length; k++) {// альфа на вектор Р і вагу ...
        				tmpVal += BaseVector.get(k)*res[k];
        			}
        			
        			for (int k = 0; k < res.length; k++) {//число на вектор Р 
        				tmp[k]+=BaseVector.get(k)*tmpVal*val;
        			}
			}
        	val = func.evaluateDerivaty((x-xii)/rho);
        	BaseVector = this.getBaseVector(xii,yii, xii,yii);//getBaseVector(xii, xii);
            for (int i = 0; i < b.length; i++) {
                b[i] = BaseVector.get(i)*val - tmp[i];
                if(Math.abs(b[i])>EPS){key = true;}
            }
            if(key){
            	double[] res2 = lu.calculate(b);
            	BaseVector = this.getBaseVectorDerivatyX(x,y, x,y);//this.getBaseVectorDerivaty(x, x);
            	for (int j = 0; j < res2.length; j++) {
					returnValue+= BaseVector.get(j)*res[j];
				}
            	BaseVector = this.getBaseVector(x,y, x,y);//this.getBaseVector(x, x);
            	for (int j = 0; j < res2.length; j++) {
					returnValue+= BaseVector.get(j)*res2[j];
				}
            }
            
            
        }
        return returnValue;
    }

    public double calculateDerivatyY(double x, double y) {
        throw new UnsupportedOperationException("Not supported yet.");
    }

    public static void main(String[] args){
        ArrayList<Particle> nodes = new ArrayList<Particle>();
        ArrayList<JMeshlessInterpolationFunction> listBaseFunction = new ArrayList<JMeshlessInterpolationFunction>();
        //параметри області:
        double ax = 0, bx = 1,
               ay = 0, by = 1;

        int n = 25, m = 25;

        double  hx = (bx-ax)/(double)n,
                hy = (by-ay)/(double)m;

        //вузли+функції:
        for (int i = 0; i <= n; i++) {
            for (int j = 0; j <= m; j++) {
                Particle p = new Particle(ax+i*hx, ay+j*hy, i*n+j);
                nodes.add(p);
            }
        }
        final JMeshlessInterpolationFunction2D ff = new JMeshlessInterpolationFunction2D(1, 0.25, 0.25, 0.5, new JGaussian(0.5));
        ff.setNodeList(nodes);
        Mapper mapper = new Mapper() {
			
			@Override
			public double f(double x, double y) {
				return ff.calculate(x, y);//(x, y);
			}
		};
		final JMeshlessInterpolationFunction2D ff1 = new JMeshlessInterpolationFunction2D(1, 0.5, 0.5, 0.5, new JGaussian(0.5));
        ff1.setNodeList(nodes);
		Mapper mapper1 = new Mapper() {
			@Override
			public double f(double x, double y) {
				return ff1.calculate(x, y);
			}
		};
		
		
		Range xrange = new Range(0.0,1.0);
		Range yrange = new Range(0.0,1.0);
		int nsteps = 50;
		
		final Shape surface = (Shape)Builder.buildOrthonormal(new OrthonormalGrid(xrange, nsteps, yrange, nsteps), mapper);
		ColorMapper myColorMapper=new ColorMapper(new ColorMapRainbow(), surface.getBounds().getZmin(), surface.getBounds().getZmax(), new Color(1,1,1,.5f)); 
		surface.setColorMapper(myColorMapper);
		surface.setFaceDisplayed(true);
		surface.setWireframeDisplayed(true);
		surface.setWireframeColor(Color.BLACK);
		
		final Shape surface1 = (Shape)Builder.buildOrthonormal(new OrthonormalGrid(xrange, nsteps, yrange, nsteps), mapper1);
		ColorMapper myColorMapper1=new ColorMapper(new ColorMapRainbow(), surface.getBounds().getZmin(), surface.getBounds().getZmax(), new Color(1,1,1,.5f)); 
		surface1.setColorMapper(myColorMapper1);
		surface1.setFaceDisplayed(true);
		surface1.setWireframeDisplayed(true);
		surface1.setWireframeColor(Color.BLACK);

		// Compute an image of the contour
//		MapperContourMeshGenerator contour = new MapperContourMeshGenerator(mapper, xrange, yrange);
		
		// Create a chart with contour axe box, and attach the contour picture
		JzyFactories.axe = new AxeFactory(){
			@Override
			public IAxe getInstance() {
				return new ContourAxeBox(box);
			}
		};

		ChartMouseController mouse = new ChartMouseController();
//		Chart chart = new Chart(Quality.Nicest,"swing");
//		Chart chart = new Chart(Quality.Advanced,"swing");
		Chart chart = new Chart(Quality.Intermediate,"swing");
		chart.addController(mouse);
//		ContourAxeBox cab = (ContourAxeBox)chart.getView().getAxe();
//		ContourMesh mesh = contour.getContourMesh(new DefaultContourColoringPolicy(myColorMapper), 400, 400, 10, 0, false);
//		cab.setContourMesh(mesh);
		
		// Add the surface and its colorbar
//		chart.addDrawable(surface);
		chart.addDrawable(surface1);
		surface.setLegend(new ColorbarLegend(surface, 
				chart.getView().getAxe().getLayout().getZTickProvider(), 
				chart.getView().getAxe().getLayout().getZTickRenderer()));
		surface.setLegendDisplayed(true); // opens a colorbar on the right part of the display
		
		
		
		JFrame frame = new JFrame();
		frame.add((javax.swing.JComponent)chart.getCanvas());
		frame.setSize(new Dimension(400, 400));
		frame.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
		frame.setVisible(true);
        
        for (int i = 0; i < nodes.size(); i++) {
            Particle particle = nodes.get(i);
            System.out.println(particle.x+"\t"+particle.y+"\t"+ff.calculateDiffuseDerivatyY(particle.x, particle.y));
        }

    }


//    public double calculateDerivaty(double x){
//        ArrayList<Double> BaseVector = getBaseVector(x, xii);
//        int n = BaseVector.size();
//        b = new double[n];
//        A = new HashMatrix2(HashMatrix2.ARRAY_TYPE, n);
//        for (int i = 0; i < NodeList.size(); i++) {
//            Particle particle = NodeList.get(i);
//            calculateAlfa(x,particle.x);
//        }
//
//        boolean key=false;
//        for (int i = 0; i < b.length; i++) {
//            b[i] = BaseVector.get(i)*func.evaluate((x-xii)/rho);
//            if(b[i]!=0.0){key = true;}
//        }
//
//        double returnValue = 0.0;
//        if(key){
//            LUFactorization lu = new LUFactorization(A);
//            double[] res = lu.calculate(b);
//
//            BaseVector = getBaseVectorDerivaty(x, xii);
//
//            for (int i = 0; i < b.length; i++) {
//                returnValue+=BaseVector.get(i)*res[i];
//            }
//        }
//        return returnValue;
//    }

}
